A thermodynamically reversible generalization of Diffusion Limited Aggregation 
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We introduce a lattice gas model of cluster growth via the diffusive aggregation of particles in 
a closed system obeying a local, deterministic, microscopically reversible dynamics. This model 
roughly corresponds to placing the irreversible Diffusion Limited Aggregation model (DLA) in con- 
tact with a heat bath. Particles release latent heat when aggregating, while singly connected cluster 
members can absorb heat and evaporate. The heat bath is initially empty, hence we observe the 
flow of entropy from the aggregating gas of particles into the heat bath, which is being populated 
by diffusing heat tokens. Before the population of the heat bath stabilizes, the cluster morphology 
(quantified by the fractal dimension) is similar to a standard DLA cluster. The cluster then grad- 
ually anneals, becoming more tenuous, until reaching configurational equilibrium when the cluster 
morphology resembles a quenched branched random polymer. As the microscopic dynamics is in- 
vertible, we can reverse the evolution, observe the inverse flow of heat and entropy, and recover the 
initial condition. This simple system provides an explicit example of how macroscopic dissipation 
and self-organization can result from an underlying microscopically reversible dynamics. 

We present a detailed description of the dynamics for the model, discuss the macroscopic limit, 
and give predictions for the equilibrium particle densities obtained in the mean held limit. Empirical 
results for the growth are then presented, including the observed equilibrium particle densities, the 
temperature of the system, the fractal dimension of the growth clusters, scaling behavior, finite size 
effects, and the approach to equilibrium. We pay particular attention to the temporal behavior 
of the growth process and show that the relaxation to the maximum entropy state is initially a 
rapid non-equilibrium process, then subsequently it is a quasistatic process with a well defined 
temperature. 



PACS number(s): 05.70.Ln, 05.70.-a, 61.43.Hv, 



MICROSCOPIC REVERSIBILITY AND 
PATTERN FORMATION 



Pattern formation is an intrinsically dissipative process 
11, however the laws of physics are microscopically re- 
versible: there is no dissipation at the microscopic scale. 
In this paper we describe a simple system which organizes 
into patterns through microscopically reversible dynam- 
ics, hence it also models how dissipation arises (i.e., how 
information flows between the macroscopic and the mi- 
croscopic degrees of freedom). This system provides a 
clear example of how to reconcile the macroscopic irre- 
versibility that gives rise to patterns with the microscopic 
reversibility adhered to by physical processes. Motivated 
by the desire to understand this general issue, we study 
specific details of the model, focusing on transitions in 
the resulting growth morphology and the approach to 
thermodynamic equilibrium. 

We have previously observed several examples of re- 
versible cellular automata dynamics which produce large 
scale order through microscopically reversible dynamics 
[§-f|]. In contrast other research in the field of pat- 
tern formation has focused on irreversible microscopic 
mechanisms, with examples ranging from crystal growth 
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||[| , to Turing patterns in chemical reactions , to pat- 
terns formed by growing bacterial colonies || , to kinetic 
growth problems |)|,|l0| . 

Here we model cluster growth by reversible aggrega- 
tion (RA), generalizing the irreversible Diffusion Lim- 
ited Aggregation model (DLA) [jlO| to include contact 
with a heat bath. Particles, which are initially diffus- 
ing on a two-dimensional lattice, stick upon first con- 
tacting a cluster member and release heat which then 
diffuses about a superimposed lattice representing the 
heat bath. The two subsystems exchange only heat and 
together form a closed system. The release of heat trans- 
fers entropy from the aggregating system (which is be- 
coming ordered) into the heat bath (which was initially 
empty). When the heat bath is nearly empty the model 
is essentially equivalent to the canonical DLA formula- 
tion (analogous to a supercooled gas crystalizing in a far 
from equilibrium situation). Hence the RA growth clus- 
ter initially resembles a typical DLA cluster. As the heat 
bath becomes populated, singly connected cluster mem- 
bers are able to absorb heat and evaporate. As the ef- 
fect of evaporation becomes significant the RA and DLA 
models diverge. The RA dynamics is exactly invertiblc: 
at any point we can invert the dynamics and run back- 
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wards, observing the flow of heat from the heat bath back 
into the gas-crystal system until we recover the exact ini- 
tial condition. 

The population levels of the heat bath and of the ag- 
gregate initially grow linearly in time, quickly reaching 
stable values which remain very nearly constant for the 
remainder of the evolution. The energy of each subsys- 
tem is a function only of the population levels, indepen- 
dent of the physical configuration of the particles. Hence, 
once the population levels stabilize, the rate of energy ex- 
change (which is entirely in the form of heat) between the 
heat bath and the gas-aggregate system becomes so slow 
that we can characterize the subsequent dynamics as a 
quasistatic process, with a well defined temperature at 
all times. 

The aggregate mostly forms while the heat bath is at 
a lower temperature than in the quasistatic steady-state. 
Hence, after the population levels stabilize, the clus- 
ter slowly anneals. The cluster morphology, which ini- 
tially resembles a DLA cluster, gradually becomes more 
spread out and tortuous, until it ultimately resembles a 
branched polymer with quenched randomness. The two 
timescales that characterize the growth process are sep- 
arated by two orders of magnitude. Initially, the popula- 
tion levels quickly reach a quasistatic steady-state. Sub- 
sequently, the aggregate slowly anneals until reaching the 
ensemble of configurations corresponding to the highest 
entropy macrostate (the branched polymer). 

Aside from insight into microscopically reversible 
mechanisms that give rise to macroscopic patterns, the 
development of invertible dynamics and algorithms has 
technological significance in pushing down the barrier to 
atomic scale computing. Each bit of information erased 
at temperature T releases at least TAS = fcsTln2 units 
of heat into the environment pl| . Heat is created in pro- 
portion to the volume of the computer, yet heat leaves the 
computer only in proportion to the surface area. Hence, 
as logic gate density in computers increases, the use of 
an invertible dynamics (which does not erase informa- 
tion and hence does not need to produce heat) will be 
required to keep the mechanical parts from burning up 
p2| , ^3| . From a more pedagogical viewpoint, discrete 
computer models of reversible microscopic dynamics pro- 
vide a laboratory for studying non-equilibrium statistical 
mechanics and the approach to equilibrium. These mod- 
els let us explore physically plausible dynamics for non- 
equilibrium systems (i.e., discrete dynamics which are 
microscopically reversible and thus automatically obey 
Liouville's Theorem). A particularly instructive example 
of this approach is the formulation of a dynamical Ising 
model gj -fi6[. However, more widely used in physics 
arc discrete, reversible models of fluid flow such as the 
HPP and FHP lattice gases 0Q. For a recent dis- 
cussion of modeling physical phenomena with reversible 
computer models see Ref. || . For a recent discussion of 
macroscopic irreversibility and microscopic reversibility 
see Ref. jF9|. For a recent discussion of techniques for 
the explicit construction of reversible models in statisti- 



cal mechanics see Ref. [£0| ; but note that closely related 
techniques were discussed in the early 1980's by Fredkin 
(as discussed in Ref. jpjl)- 

The initial sections of this manuscript describe our 
model; the middle contain a mathematical formulation 
of the model; the final, the empirical results. Specifi- 
cally, Sec. H describes the detailed dynamics, includ- 
ing the subtleties of constructing an invertible dynam- 



ical model and implementation issues. In Sec. Ill 



we 

discuss the macroscopic limit of an analytic formulation 
of the model and establish the reaction-diffusion equa- 
tions describing the system. In Sec. [V we treat the 
reaction-diffusion equations in the mean field limit and 
compare predictions for equilibrium densities of particles 
to empirical measurements. Empirical measurements of 
temperature are presented in Sec. VA, with emphasis 



on the quasistatic nature of th e annealing portion of the 
growth process. In Sec. VB we study the evolution of 



the fractal dimension of the clusters and thus quantify 
the change in growth morphology as the clusters relax 
to the maximum entropy state. We conclude with a dis- 
cussion of limitations and possible modifications of our 
model. 



II. MODELING AGGREGATION 

A. Diffusion Limited Aggregation 

Diffusion Limited Aggregation (DLA) Q is a concep- 
tually simple model which serves as a paradigm for some 
aspects of kinetic growth phenomena. Several compre- 
hensive reviews of DLA have been written. In particular 
see Ref. [^2) for a clear presentation of the basics, Ref. |2j| 
for details on physical mechanisms, and Ref. [£4| for a re- 
view of real-space renormalization group approaches to 
DLA. 

The typical scenario for DLA begins with a vacant two- 
dimensional lattice initialized with a single stationary 
seed particle, which is the nucleation site for a growth 
cluster. Moving particles are introduced from the edges 
of the lattice, following random walks along the lattice 
sites. When a moving particle lands on a site adjacent to 
a stationary seed particle (an active site) it sticks (i.e., 
the moving particle aggregates and becomes a stationary 
seed particle). The frozen aggregate particles constitute 
the solid (crystal) phase, and the moving particles con- 
stitute the gas phase. Aggregation hence consists of a 
particle undergoing an irreversible transition from gas 
to solid. Gas particles are usually introduced in a se- 
rial manner: only one gas particle is diffusing at a time. 
However, to take advantage of parallel computational re- 
sources, parallel models of DLA have been studied in 
which multiple particles are diffusing at once p5| , p6[ . In 
the dilute particle limit these models recover the serial 
DLA model exactly. 
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With the first aggregation event the DLA cluster grows 
from one single to two adjacent sites. The presence of 
the second cluster member eliminates certain paths along 
which random walkers could approach the first, with the 
effect that the probability of sticking at either end (tip) of 
the cluster is enhanced, whereas the probability of stick- 
ing along the edge of the cluster is reduced. As particles 
continue to aggregate creating new cluster tips and edges, 
the probability to stick at the tips continually outweighs 
the probability to stick along the edges. This leads to 
branching. A second influence on the growth morphol- 
ogy comes from shadowing: the probability for a particle 
to diffuse into the center of the growth cluster before en- 
countering an active site becomes negligible as the clus- 
ter grows in size. Hence the outer tips grow most rapidly. 
As a result the growth aggregate rapidly assumes a bushy 
and branching, random fractal structure, resembling frost 
on a window pane, the branching of neurons, and many 
other branched structures found in nature. 



B. Reversible Aggregation 

Our goal is to introduce a reversible, deterministic 
model of growth by aggregation, where reversible means 
that from any state of our system we can recover the 
previous state exactly. We must address the subtleties of 
making each component step invertible, including steps 
which realize stochastic processes. As we discuss below, 
the same mechanisms that are employed in our model 
in order to ensure exact conservation of energy, parti- 
cle number, and other constraints, also make it easy to 
incorporate invertibility. The stochastic component of 
the model is diffusion, which is modeled as a sequence of 
invertible "random walks" based on a deterministic algo- 
rithm using an invertible pseudorandom number genera- 
tor. 

Information about exactly when and where a particle 
undergoes a phase transition is stored in the heat par- 
ticles. The idea of storin g in formation in a heat bath 
was introduced by Creutz |14| to explore the connection 
between the microcanonical and canonical ensembles for 
the dynamical Ising model. Heat bath techniques have 
been used on occasion since then to construct reversible 
computer models of physical phenomena . 

1. Overview of the model 

To construct the reversible aggregation (RA) model 
we begin with a parallel DLA system (similar to that de- 
scribed above) and add degrees of freedom at each site, 
corresponding to a distributed heat bath. The latent 
heat released during each aggregation event can then be 
explicitly represented. In the RA dynamics, whenever a 
random walking gas particle lands on a site with exactly 
one nearest neighbor crystal particle, it will stick only if 



there is room locally in the heat bath to accept the la- 
tent heat it will release as the particle transitions from 
the gas phase to the crystal phase. The heat is released 
in quantized units called heat particles, with one heat 
particle released for each aggregation event. These heat 
particles diffuse amongst themselves (i.e., they undergo 
random walks along the lattice sites, independently of 
the gas particles) . Explicitly modeling the latent heat re- 
leased upon aggregating provides a mechanism for mod- 
eling the inverse process: a diffusing heat particle which 
contacts a susceptible crystal particle (a crystal particle 
which has only one nearest neighbor crystal particle) is 
absorbed while the crystal particle evaporates to become 
a gas particle which then diffuses away. 

The restriction on the dynamics that aggregation and 
evaporation events can occur only when exactly one near- 
est neighbor is a crystal particle means that only one 
crystal bond is ever formed or broken when a single lat- 
tice site is updated. As each aggregate particle con- 
tributes one crystal bond to the aggregate, and there is no 
further potential energy contribution, the energy of the 
aggregate is a function only of the number of aggregate 
particles, independent of their configuration. Moreover 
this constraint has two direct implications for the growth 
morphology. The first is that evaporation can only oc- 
cur for particles which are singly connected to the growth 
cluster, and so the aggregate cannot break off into discon- 
nected clusters. The second is that it introduces an ex- 
cluded volume (i.e., no closed loops can be formed), thus 
we might expect the equilibrium cluster configuration to 
be similar to that of a polymer. Note that introduction 
of an evaporation mechanism in the RA model mitigates 
the shadowing effect that was important in determining 
the DLA growth morphology: crystal particles within the 
cluster can evaporate, thus introducing gas particles into 
the interior of the aggregate. 



2. The detailed dynamics 

The RA model is constructed with 7-bits of state at 
each lattice site. One bit, N c (x,t), denotes the pres- 
ence or absence of a crystal particle at that site (i.e., 
N c (x,t) = 1 indicates presence, N c (x,t) — absence). 
Two bits, N^(x,t) where 7 = {1,2}, denote the pres- 
ence or absence of each of two gas particles. Two bits, 
N7(x,t), denote the presence or absence of each of two 
heat particles. The final 2 bits, £ s (c?, t), £h(%, t), are inde- 
pendent binary pseudorandom variables. The dynamics 
of the model consists of two kinds of steps: diffusion steps 
alternating with interaction steps. 

The same kind of diffusion process is applied to the gas 
and heat subsystems simultaneously and independently, 
while the crystal particles remain unchanged. A given 
diffusion step consists of two parts: mixing and trans- 
port. During the mixing portion of the step, a binary 
random variable is used to determine whether or not the 
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two particle bits of that species at the site (x, t) are in- 
terchanged: 

N? = (l-£ i )N? + t i Nt, (1) 

where i = g or i = h. During the transport portion of 
the step, every site replaces its first particle bit (7=1) 
with the first particle of its neighbor a distance k away 
on one side, and its second particle (7 = 2) with the one 
from the same distance neighbor on the opposite side. 
At even time steps, we use horizontal neighbors (i.e., the 
diffusion moves particles horizontally): 

Nl(x,t + 1) = N?(x + kx,t) 

N?(x,t+l) = N?(x-kx,t). (2) 

At odd time steps we use vertical neighbors (i.e., substi- 
tute y for x in Eq. (§)). The only differences between 
the gas and heat diffusion are (1) each uses a separate 
binary random variable to control its mixing, and (2) the 
distance of the neighbor particle to be copied, k, can be 
chosen separately for each subsystem — this allows us to 
independently control the diffusion constants for the heat 
bath and for the gas (c.f. Q). 

Diffusion steps alternate with steps in which the two 
subsystems — gas-crystal and heat bath — interact allow- 
ing aggregation and evaporation. The rule at a single 
lattice site during an interaction step is that exactly one 
particle can aggregate or evaporate provided that 

(a) there is exactly one crystal particle at one of the four 

nearest neighbor sites 

(b) there is room at the site for a crystal particle (for 

aggregation) or for another gas particle (for evapo- 
ration), and 

(c) the heat bits at the site can absorb (for aggregation) 

or supply (for evaporation) one unit of heat. 

Since the gas and heat particles will undergo a mixing 
step before transport, it makes no difference which of the 
two available gas particle positions a crystal particle is 
moved into when it evaporates, or which of the two pos- 
sible heat particle positions a unit of heat gets put into. 
Defining this precisely will, however, become important 
when we discuss invertibility. 

The interaction rule described thus far would be suf- 
ficient if we updated just one lattice site at a time. If, 
however, all sites on the lattice are updated simultane- 
ously, then the global dynamics no longer obeys the "sin- 
gle bond" constraint — that at any site where particles ag- 
gregate or evaporate exactly one crystal bond is formed 
or broken. For example, suppose that the tip of a crystal 
branch evaporates just as a gas particle condenses next 
to it. Each of these events would separately preserve 
the constraint, but the two simultaneous events result 
in the addition of a disconnected crystal particle which 



has no other crystal particle immediately adjacent to it. 
We can easily avoid this difficulty by holding the values 
at the adjacent sites fixed during a step in which we let 
the subsystems interact at a given lattice site, since the 
interaction step has a nearest neighbor range. In other 
words, we perform a checkerboard updating: all of the lat- 
tice sites in which the x and y lattice coordinates add up 
to an even number are updated using our single site in- 
teraction rule, while the odd sites are held fixed, and then 
vice versa. Since nearest neighbors are held fixed during 
an interaction, the constraint is obeyed locally and thus 
it is also obeyed globally. The overall dynamical rule is 
summarized in Table | — the various phases of the rule 
are applied consecutively. 

1. interact gas/heat/crystal at even sites 

2 . mix gas and mix heat separately 

3. transport gas and heat horizontally 

4. interact gas/heat/crystal at odd sites 

5 . mix gas and mix heat separately 

6. transport gas and heat vertically 

TABLE I. The various phases of one step of the RA dy- 
namics. Each phase is applied over the entire lattice simulta- 
neously. 
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Every phase of the rule described in Table [J can be 
inverted. The transport portion of the step can be run 
backwards by simply moving all particles back into the 
sites they came from (i.e., inverting the directions of the 
transport in Eq. |^). The mixing portion of the step is 
easy to invert, given the same "random" binary vari- 
ables that were used to determine which pairs of bits 
were originally swapped. We simply swap exactly those 
pairs again. The pseudorandom portion of the system 
(which supplies the random data) can simply follow some 
invertible dynamics that is independent of the rest of the 
system — the rest of the system looks at the state of this 
subsystem but does not affect it — so this pseudorandom 
subsystem can be run backwards independently of every- 
thing else. 

Making the interaction steps invertible is also straight- 
forward. When a single gas particle turns into a crystal 
particle, we put the heat token that is released into the 
corresponding heat bit (i.e., the heat particle with the 
same value of 7), and thus we remember which of the 
two gas particle bits was initially occupied. If the corre- 
sponding heat bit is already occupied, the particle is not 
allowed to aggregate (even if the other heat bit is unoccu- 
pied) . Similarly, a crystal particle is allowed to evaporate 
only if it can move into the gas bit with the same value 
of 7 as the heat token being absorbed. If there are two 
gas particles at a site we impose the constraint that the 
particle with 7=1 attempts to aggregate first; likewise if 
there are two heat particles at a site, the crystal particle 
attempts to absorb the 7=1 particle first. This does 
not introduce a bias to the growth since we are mixing 
the 7=1 and 7 = 2 variables in an unbiased manner at 
each time step. With these refinements, our interaction 
rule applied to a single site is its own inverse: if we apply 
it twice (without a diffusion step in between) we get back 
the state we started with. Since the interaction rule is 
applied in a checkerboard fashion, sites are updated in- 
dependently: if we apply the rule a second time to the 
same checkerboard, it will undo the first application at 
every site. 

Thus an inverse step consists of applying the inverses 
of the rule-phases in the opposite of the order listed — 
once one phase is undone, the previous phase can be 
undone. Each inverse step undoes one step of the for- 
ward time evolution. As we watch the inverse evolution, 
we see each heat particle retrace its path, to be in ex- 
actly the right location at the right time to uncrystalize 
the crystal particle which originally released it. Particles 
un-aggregate and un-diffuse and un-evaporate in a man- 
ner that exactly retraces their behavior in the forward 
evolution. 



3. Implementation 

The RA model was implemented on a special pur- 
pose cellular automata machine, the CAM-8 ^9|, which 



was designed to efficiently perform large-scale uniform, 
spatially arrayed computations. On this machine, all 
simulations must be embedded into a lattice gas frame- 
work [^,^|, in which uniform data movement (data- 
advection) alternates with processing each site indepen- 
dently (site- update). For a 2-dimensional model such as 
ours, sheets of bits move coherently during the advec- 
tion phase: corresponding bits at each site all move in 
the same direction by the same amount. The boundaries 
are periodic — bits that shift past the edge of the lattice 
reappear at the opposite edge. After moving the bits, we 
perform the site update phase. During this phase, the 
bits that have landed at each lattice site are updated in 
a single operation by table lookup: the bits at each lat- 
tice site are used as an index into a table that contains a 
complete listing of which state should replace each pos- 
sible original state. Both the data movement and the 
lookup table can be freely changed between one lattice 
updating step and the next. 

Our model requires 7 bits of state to appear at each 
site in our LxL lattice. Using random data generated by 
a serial computer, the bits which correspond to the gas 
particles are initially randomly filled with a 4% density 
of particles and the bits which correspond to the binary 
random variables with a 50% density of particles. One 
crystal particle is placed at the center site of the lattice. 
The heat bath is initially empty. 

The dynamics on the pseudorandom subsystem is very 
simple: each of the two random bit-planes (each consist- 
ing of all the ?y g 's or all of the r/hs) are simply shifted by 
some large amount at each time step. We could choose 
the amount and direction of each shift at random for 
each step of updating, using a reversible random num- 
ber generator running on the workstation controlling the 
simulation. Instead, the simulations discussed here sim- 
ply shift the bit planes by a large and fixed amount at 
each step, making sure that the the x and y components 
of the two shifts are all mutually co-prime, as well as be- 
ing co-prime with the overall dimensions of the lattice. 
Thus to run the random subsystem backwards, we just 
reverse the direction of the shifts. 

The checkerboard updating is accommodated by 
adding an eighth bit to each lattice site, and filling these 
bits with a checkerboard pattern of ones and zeros. In our 
rule the various subsystems are allowed to interact only at 
sites marked with a one. To change which checkerboard 
is marked for updating, we simply shift the checkerboard- 
marker bit-plane by one position in the +x direction. 

The rule described in Table Q turns into two lattice-gas 
steps on CAM-8. The first three phases listed in Table || 
arc done during one step, and the next three in the second 
step. The data movement is part of each step: note that 
each of the two steps uses the same lookup table applied 
to each lattice site, but slightly different data movement. 
To run backwards, we use the inverse lookup table, and 
the inverse data movement. Note that in the discussion 
of experimental results, everything in Table || is counted 
as a single step. 
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CAM-8's event counting hardware was used to mon- 
itor simulation parameters while the simulations ran. 
Including event counting, the 8-processor CAM-8 per- 
formed about 10 8 site update operations per second for 
this model. 



III. THE MACROSCOPIC LIMIT 

The dynamics of the RA model, described in detail 
above, can be succinctly presented in an analytic frame- 
work. We develop this framework first in terms of the 
discrete space, time, and occupation number variables. 
We then ensemble average over the occupation numbers 
and take the continuum limit of the space and time vari- 
ables to establish the reaction-diffusion equations for the 
system. 

As discussed in Sec. II B 2| there are seven bits of state: 
N c (x,t),N](x,t),N^(x,t),£g(x,t), and ^n(x,t), where 
7 = {1,2}. They correspond respectively to one bit of 
crystal, two bits of gas particles, two bits of heat par- 
ticles, and two bits of random data. N?(x,t) = 1 in- 
dicates the presence of species i at location x and time 
t, in channel 7, and N?(x,t) = indicates the absence. 
The absence or presence of a crystal particle is denoted 
by N c (x,t) — {0,1} respectively. The total number of 
particles of species i at time t present on the lattice is 
denoted by Mi{t) = t)i where the sum is over 

all of the lattice sites and the two particle channels. 

There are no external sources or sinks for any of the 
three species represented (the gas, crystal, and heat 



species). Energy is only exchanged between the gas- 
crystal and the heat bath subsystems. Thus the complete 
system is isolated. Conservation of the total number of 
gas and crystal particles implies that 



Af g (t)+M c (t)=M g (0)+N c (0)- 



(3) 



Conservation of the total energy of the system implies 
that 



N g {t)e g +N c {t)e c 



M h {t)e h -- 



-K{0)e c +M h (0)e h , (4) 



where £j represents the energy (kinetic and potential) 
per particle of species i (notice that there is no config- 
urational contribution to the energy of the crystal). As 
discussed in Sec. |0| each aggregation event releases one 
heat particle (likewise each evaporation event absorbs one 
heat particle), thus Eh — s g ~ £ c and moreover 



(5) 



(Note that Ah(0) = and A/" c (0) = 1 in our experiments). 

To facilitate the description of the dynamics, we intro- 
duce a functional equation which is +1 at any site where 
a particle is about to crystalize, —1 at a site where a par- 
ticle is about to evaporate, and otherwise. The func- 
tional, T~< [N](x,t),iq(x,t),N c (x,t),{N c (x + e k ,t)}], 
is evaluated on a neighborhood of lattice sites surround- 
ing some given position x at a given time t (the notation 
{N c (x + Bk,t)} refers to the set of values of N c for the 
nearest neighbors of the point x). 



T^(x, t) = N]{x, t)[l - Nl{x, i)][l - N c {x, t)] Y / N c (x + e h t) J] [1 - N c (x + e k ,t)\ 

3=1 kjtj 
d 

-[1 - N](x, t)]N%(x, t)N c (x, t) J2 N c (x + ej,t) N c (x + e k , t)} . (6) 



Here ij and e k are the vector lattice directions of the 
nearest neighbor cells, and d the number of distinct lat- 
tice directions. For a two dimensional square lattice (i.e., 
the lattice used for the present implementation) d = 4 
and the vector lattice directions are {£, —x,y, —y}- 

The first term in Eq. (|^) equals 1 if a gas particle in 
channel 7 is present at site x and time t, a heat par- 
ticle in channel 7 is absent at site x and time t, there 
is no crystal particle already at that site, and only one 
crystal particle is present at a nearest neighbor site. It 
is zero otherwise. The second term equals 1 if there is 
no gas particle in channel 7 present at site cc and time 
t, there is a heat particle in channel 7 present at site 
x and time t, there is a crystal particle present at that 
site, and only one crystal particle is present at a nearest 



neighbor site. It is zero otherwise. The first and second 
terms are mutually exclusive (a heat particle in channel 
7 cannot be simultaneously present and absent, nor can 
a gas particle). 

The dynamics consists of making the changes indicated 
by T 1 and then T 2 , then applying a random permutation 
to mix 7=1 and 7 = 2, and then performing the stream- 
ing step to move the particles. The permutation mixes 
the N} (x, t) and A^ 2 (x, t) components in an unbiased way, 
so it is simpler to discuss the dynamics of a combined 
variable, Ni(x,t) — Nl(x,t) + Nf(x,t). Likewise, if we 
let T{x, t) = [T 1 (x , t) + T 2 (x , t)\ [1 - T 1 (£, t)T 2 (x , t) /2] , 
the interaction portion of the dynamics at a single lattice 
site can be written 

N c (x, t + 1) = N c (x, t) + F{x, t) 



G 



Ng{x, t+l) = Ng{x, t) ~ t) 

N h (x, t + l) = N h (x, t) + T{x, t) 



(7) 



The [1 - J 7l {x,t)J 72 (x,t)/2] factor in the definition of 
T(x,t) reflects the fact that only one transition occurs 
at a given site, even if two gas or two heat particles are 
present. 

To construct continuous variables from the discrete 
ones, we consider the average short-term behavior of 
the system over an ensemble of independent realizations 
which all have the same set of local particle densities. For 
each discrete variable, we let rii(x, t) = (Ni(x, t)), and de- 
note the average of the functional as (T) . (This technique 
of averaging over many independent realizations, i.e., es- 
tablishing the one particle density function, is commonly 
used to derive the lattice Boltzmann equation starting 
from discrete particle models of hydrodynamics pO| , pT[ |). 

With this notation, the average propagation of the gas 
and heat particles can be expressed as 

rii(x,t+ 1) = 



- [n,-(a? — x,t) + m(x - y, t)+ 
+m(x + x, t) + m(x + y, t)] . 



(8) 



To establish the continuum limit, we Taylor expand. The 
terms involving the first derivatives cancel, leading to the 
result 

rii(x, t + 1) = rii(x, t) + 

E^Sm^)) + o(Ax 3 ) 



4 dx 2 



rii(x,t) + 



-V 2 ni(f, t), 



(9) 



where i = g or i = h. Note that to order At the above 
equation is the standard diffusion equation, 



-grn{x,t) = ^^V 2 n 4 (x,t). 



(10) 



As discussed in Sec. IIB2, we can control the length 
of each diffusion step separately for the heat and for the 
gas particles. The heat particles execute random walks 
composed of individual steps of length k, whereas the gas 
particles execute walks of step length unity. Thus if the 
| Ax | 2 that appears in Eq. M) and in Eq. (|l^) refers to the 
gas subsystem, then \Ax\ for the heat subsystem (and 
hence its diffusion constant) is a factor of k 2 larger. 

To proceed further, we will make the approximation 
that our average variables are independent. Then we 
are allowed to replace the average of a product of vari- 
ables by the product of the average for each variable: 
(ab) = (a)(6), for a, 6 independent variables. This is 
the assumption of molecular chaos, which is also used 
to derive the lattice Boltzmann equation. With this ap- 
proximation the average of the functional J-" 1 is simply 



= {nj(x, t) [1 - nl(x, t)] [1 - n c (x, t)] 
- I 1 -rig&t)] nl[x,t)n c (x,t)} 
x E 3 - n c{x + ej, t) U k7tj [1 - n c (x + e k ,t)] ■ (11) 

Similarly, we can write down an expression for (J-) . 

To obtain the continuum limit of these averaged equa- 
tions, we again use Taylor series approximations. In the 
diffusive regime, At ~ (A/) 2 , so we truncate the expan- 
sions at these appropriate orders. Let T be the contin- 
uum limit of (J-) (which we won't write out explicitly). 
Then from Eq. (0), we obtain 



9 ,~ ^ 1 * 

-n c {x,t) = -F, 



(12) 



The other reaction-diffusion equations for our system 
can be obtained by proceeding as we did in Eq. (|9|). For 
example, under the full dynamics (which consists of both 
the interaction and diffusion phases) 



n g (x, t + l) + n c (x, t + 1) = n c (x, t) + 
\ [n g (x -x,t) + n g (x - y, t)+ 
n g (x + x,t) + n g (x + y,t)] , 



(13) 



since any particles present at a site at time t + 1 were ei- 
ther already there at time t, or moved there. Expanding 
this exactly as in Eq. (||) and using Eq. (|lj), we get 



g^n g {x,t) 
and similarly, 



d_ 

St 



nh{x, t) 



\Ax\ 2 „ _ . d _ , 
^-V 2 n g {x 1 t)--n c (x,t) 



k 2 \Ax\ 2 ^ 2 _ . d . 
4At V 2 n h (x,t) + -n c (x,t). 



(14) 



(15) 



Note that if k = 1, and we add the last two equa- 
tions the variable n g {x,t) + n,h(x,t) obeys the diffusion 
equation, unaffected by the interaction between the sub- 
systems (i.e. if we remove the distinctions between gas 
and heat, the combined variable simply diffuses without 
interacting). 

To test the consistency between the microscopic dif- 
fusive dynamics implemented in our model and macro- 
scopic descriptions given byEq. (|TJ) and Eq. @, we em- 
pirically measured the diffusion coefficient for gas and for 
heat particles as they diffuse about the space. Each par- 
ticle should execute a random walk. The variance of the 
distance from the origin in the x or y direction, a 2 , is pro- 
portional to the diffusion coefficient in that direction, Di, 
where i = g or i = h. The exact relation is Di = af/4p, 
where p is the number of steps taken. For an unbiased 
random walk, the variance of the net displacement from 
the origin is of = k 2 p, thus Di = k 2 /A. For the gas par- 
ticles (k = 1) we find D g = (0.996 ± 0.009)/4. For the 
heat particles, with k = 3, we find D h = (9.00 ± 0.08)/4. 
Thus the ratio of the heat to the gas diffusion length 
Dh/D g — 9.0 ± 0.1, agreeing with the theoretically pre- 
dicted value of k 2 . 
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IV. THE MEAN FIELD LIMIT 



The mean field limit corresponds to the "well-stirred 
reaction," meaning that the reacting species are uni- 
formly spread throughout the space, and thus each par- 
ticle feels the presence of the mean concentration of each 
species. For our system in equilibrium the gas particles 
and the heat particles are uniformly distributed through- 
out the space; it is only the crystal particles which do not 
obey this assumption. A uniform distribution means that 
there are no concentration gradients (Vn,i(x, t) = for all 
x and i), and thus X7 2 rii(x, t) — 0. Also we can drop the 
explicit x notation from the argument of the variables: 
Ui(x,t) = ni(t). 

Once the population level of the heat bath has reached 
the quasistatic steady-state, the concentrations of all 
three species will remain essentially constant and the sys- 
tems will have a well defined temperature from then on, 
as discussed in Sec. V A. We denote the time to reach the 
steady-state (i.e. the time for the subsystems to reach 
the same temperature) as tt- We can now drop the 
explicit time notation from the arguments of the vari- 
ables in steady-state: rii(t > tt) = rij. In this regime 
drii/dt = and likewise dnj/dt = 0, thus < T 1 >= 0: 



= (J^) 

= 4 [rig (1 - n h ) (1 - n c ) - (1 - n g ) n h n c ] x 
n c (1 - n c y 



OL=128 
□ L = 256 
©L = 512 



0.1 
N10) I L 2 



1. The empirically determined equilibrium value of 

2 



FIG 

n c /(l — n c ) as a function of the initial gas density, Af g (0) / L 
for systems of size L = 128, 256, and 512. The solid line is 
the mean field prediction. Note that the error bars are the 
same size as the points. 



The mean field approach makes predictions about the 
overall density of the system (hence the equilibrium tem- 
perature, as described below), but it does not make pre- 
dictions about the growth morphology. 



V. EMPIRICAL ANALYSIS 



(16) 



A. Temperature 



There are three solutions to Eq. (16). Each solution 



corresponds to fixed point of the dynamics. Only one 
is in the regime of interest. The fixed point at n c — 
corresponds to the presence of only gas particles. The 
fixed point at n c = 1 is not allowed by the aggregation 
conditions (the aggregate can not have any closed loops). 
The remaining fixed point predicts that the equilibrium 
condition is n g [1 — n^] [1 — n c ] = [1 — n g ] riyjic. Noting 
that in the mean field limit ensemble averages equal spa- 
tial averages (i.e., n g (t) = Nn(t)/ L 2 ), the constraints de- 
scribed in Eq. (|J) and Eq. (g) can be written respectively 
as n g (t) = JV g (0)/L 2 - n c (t) + l/L 2 « JV g (0)/L 2 - n c [t), 
andn c (£) = rih(t) + l/L 2 ss n/,. After incorporating these 
relations the equilibrium condition can be expressed as 



1 - n c 



■■M g (oy 



0{nl). 



(17) 



Figure [l] is a plot of the equilibrium value n c /(l — n c ) 
versus the initial density of gas particles, J^ a (0) /L 2 . The 
solid line is the mean field prediction, Eq. (|17|). The 
points were obtained empirically from our simulations 
of three different system sizes, L = 128, 256, and 512. 
The agreement between the three system sizes should be 
noted. 



The gas- aggregate subsystem and the heat bath sub- 
system together form a thermodynamically isolated sys- 
tem. These two subsystems are allowed to exchange en- 
ergy only between themselves, and this energ y is pu rely 
in the form of heat (AQ). As discussed in Sec. II B 1 and 
Sec. [II the energy of the aggregate is a function only of 



the number of aggregate particles and is independent of 
their configuration. The total internal energy of the gas 
and heat particle species also is a function only of the 
number of particles of each species. Hence if invariant 
average population densities are achieved there is no fur- 
ther net exchange of heat between the subsystems, and 
they have by definition attained the same temperature. 

The standard expression for the temperature of a two- 
level system |52| , such as the heat bath in the RA model, 
follows directly from combining the definition of temper- 
ature (1/T = AS/ AE\ V ) with the microcanonical defi- 
nition of entropy (AS — ks ln(£!//17i), where denotes 
the number of microstates consistent with the macro- 
scopic variables): 



— - -— In ( n ' 1 
T e h \l-n h 



(18) 



Directly computing the temperature of the other subsys- 
tem is more difficult, but we can infer its temperature 
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from that of the heat bath (note the gas particles are 
free to diffuse over the crystal, hence there is no change in 
the accessible volume for the heat particles as the crystal 
changes size and conformation: the crystal does no work 
on the gas, PdV = 0). 

The approach to temperature equilibrium and a 
closeup of the subsequent fluctuations in temperature 
are shown in Fig. Figure ^(a) plots the mean oc- 
cupancy of the heat bath versus the time step into the 
simulation, with the corresponding temperature (in units 
of hsT/eh) displayed on the right vertical axis. The 
initial growth is linear, with a slope of about 1.8. It 
then levels off near the quasistatic steady-state density 
of rih — 0.031 (indicated in the figure as the dashed 
horizontal line). The results are averages over several 
independent realizations for three different system sizes, 
L = 128, 256, and 512. Note that the three systems 
reach the same steady-state densities and hence the same 
temperature, but that the time to equilibrate depends on 
the system size. The data for the three systems was col- 
lapsed onto one curve by rescaling the time axis by the 
factor L z , with z = 1.8. The time to reach the equilib- 
rium temperature is tt ~ 10L Z . Note that this scaling 
behavior has an exponent which is slightly smaller than 
the diffusion exponent: the diffu sion time is proportional 
to L 2 . As discussed in Sec. VB the fractal dimension at 
time tt ~ 1.8: the time seems to scale with the fractal 
dimensionality instead of the Euclidean dimensionality of 
the space. 

To study the details of the subsequent fluctuations we 
focus on the largest system size, 512 x 512. As mentioned, 
during the initial period the growth of the heat bath pop- 



ulation (and the size of the aggregate) is linear. It levels 
off at about 8270 particles on average, in a time which 
is less than 10 6 steps. The population continues to grow 
extremely slowly after this point, rising by an average 
of 1.8 ± 0.1 particles every 1 x 10 7 steps, as determined 
by a linear regression based on about 8000 data points 
taken at equally spaced intervals in the regime where 
t > tt- The probability that the slope is actually zero 
is 3 x 10 -8 , as determined by a t-statistic comparing the 
ratio of the obtained slope to the sum of squares differ- 
ences. Figure |](b) shows a scatter plot of every third 
of the 8000 data points, overlayed by a straight line in- 
dicating the results of the linear regression on all 8000 
points (only every third point is shown for visualization 
purposes: showing all the points results in a dense black 
cloud). The actual number of particles in the heat bath 
is indicated on the left axis, the corresponding tempera- 
ture is given on the right. Although the temperature of 
the heat bath is not constant, it is very nearly so. Once 
the population levels stabilize the subsequent dynamics 
(i.e., the relaxation to the maximum entropy state for 
the crystal) is clearly a quasistatic process. The crystal 
does continue to exchange heat with the heat bath when 
it anneals, but the net heat exchange is essentially zero 
(the net heat exchange rate is ~ 2 x 10 -7 particles per 
update of the space). 

We measured also the fluctuations in population levels 
of the heat bath for the L — 256 and 128 systems, in 
the corresponding regimes. Consistent with the L = 512 
system, we find the population rises by an average of 
2.1 ± 0.1 particles every 1 x 10 7 steps. 



0.031 

N h (t/L z )/L 2 
0.01 



10 

t /L z , z=1 .1 



L= 128 
■ L = 256 
L = 512 



- 0.25 



- 0.15 
1000 



8350 - 



8300 
N h (t) 



L= 512 

y = 8269+(1.8x10~ 7 )x 



8150 
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FIG. 2. (a) The mean density of heat bath particles as a function of time into the simulation, plotted for 
L — 128, 256, and 512. The corresponding temperature is given on the right vertical axis. The initial growth of the heat bath 
density is linear, with a slope of about 1.8. Note that the steady-state density of the heat bath, and hence the temperature, is 
equal for all three system sizes, yet the time to equilibrate scales with the system size as tt ~ 10L 1 ' 8 . (b) The actual average 
values of the total population of the heat bath as a function of time, for every third time measured beyond tt. The dashed line 
is the result of a linear least squares regression on all of the data. Note the slight drift upward with time, of about 2 particles 
per 10 7 steps. 
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B. Fractal Dimension 

The aggregate formed primarily while the heat bath 
contained less energy than its equilibrium level. Hence, 
if we continue running the dynamics, the cluster anneals; 
it evolves from a DLA-like cluster to a quenched branched 
random polymer structure. To quantify the cluster struc- 
ture we calculate the fractal dimension of the aggregate 
and specifically how the fractal dimension changes as a 
function of the time into the simulation. 

We measure the fractal dimension using a box-counting 
procedure which requires that we first establish the cen- 
ter of mass of the growth aggregate (which is typically 
not the initial seed particle, as the center of mass diffuses 
about the space as the aggregate anneals). An imaginary 
window box of edge length I is defined and centered on 
the center of mass. The number of lattice sites within 
that window that contain a crystal particle, jV c (Z), is tal- 
lied. The window size is increased and the count retai- 
ned. This procedure is iterated until the number of crys- 
tal particles contained no longer increases with window 
size. Before saturation, the number of particles contained 
should increase with some power of the window size 

Af c (l)ocl d f. (19) 

The exponent df is the fractal dimension. 

The RA cluster should initially resemble a parallel, ir- 
reversible DLA cluster of equivalent mass. Figure |^(a) 
shows a typical RA cluster for the L = 512 system at the 
time t = tt, which is the time for the mass of the RA clus- 
ter to stabilize at essentially its final mass (Af c — 8270). 
Figure ||(b) shows a typical DLA cluster of equivalent 



mass. Both systems were initialized with a 4% density 
of diffusing gas particles. The gas particles still present 
at this stage of the evolution are shown as the small dots 
in the figure. Note that for the RA system the gas par- 
ticles are distributed throughout the space, yet for the 
DLA system very few gas particles penetrate the region 
defined by the edges of the cluster. The RA cluster has 
experienced enough annealing by the time t — tt to have 
a fractal dimension different than that of the DLA clus- 
ter, yet the radii of both clusters are comparable and 
approximately equal to a quarter of the length of the 
system (r k, L/A). The RA cluster morphology is still 
far from its final morphology. 

Figure ^ shows the box-counting results obtained for 
both models in the regime described above and pictured 
in Fig. ||. The top curve is for DLA, the bottom for 
RA. Both models were implemented on a L — 512 size 
system. The vertical axis is the mass contained in the 
window, Af c (l), the horizontal is the window size I. The 
curve for the DLA system is the average of 10 indepen- 
dent DLA clusters of mass Af c ~ 8270. The curve for the 
RA system is the average of ten independent RA clusters 
sampled at time t = tt. The slope of the curve corre- 
sponds to the fractal dimension and was determined via 
a linear least squares fit. Consistent with past numerical 
studies of DLA |^2|,^3f , we find that the fractal dimension 
for DLA clusters isd° LA = 1.71 ±0.01 (for difficulties as- 
sociated with determining the fractal dimension of DLA 
see the detailed discussion in Ref . |24| ) . For the RA clus- 
ters the fractal dimension is d^ A (t = tt) = 1.81 ± 0.03. 
A line with this slope is shown overlaying each respective 
curve HJ]. 



FIG. 3. Two growth clusters of the same mass, M c ~ 8270. (a) A cluster grown via the RA model, pictured at time t = tt, 
where tt is the time for the heat bath and gas-aggregate system to reach the same temperature, (b) A parallel DLA cluster. 
Note the gas particles, which are shown as the small dots. For the RA system the gas particles are distributed throughout the 
space, yet for the DLA system very few gas particles penetrate the region defined by the edges of the cluster. 
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FIG. 4. The number of aggregate particles contained in a 
box of length I, as a function of I. The slope of the line is the 
fractal dimension. The top curve is for parallel DLA clusters 
of mass JV C — 8270. The bottom curve is for the RA clusters 
sampled at t = tj?. Examples of these clusters are pictured in 
Fig. | 



The RA cluster is less dense than the DLA cluster in 
the area immediately surrounding the initial aggregation 
site, however the radii of both clusters are comparable. 
Several of the aggregate particles in the RA cluster have 
annealed away from the center to occupy the region be- 
tween the center and the edge of the cluster. Hence, at 
the time depicted in Fig. ||, the RA cluster has a higher 
fractal dimension than the DLA cluster. The only con- 
straints on the cluster are the number of particles and 
the connectivity. As there are more ways to have con- 
nected clusters of a specified particle number in an area 
of larger radius, the RA cluster evolves from the dense, 
bushy DLA-like structure shown in Fig. to a tenu- 
ous structure which occupies more of the available lat- 
tice (with the initial increase in fractal dimension being 
a transient behavior). We ultimately expect to observe a 
diffuse structure with just a few meandering vines which 
can access more of the available configuration space. 

As the time into the simulation advances, the density of 
the growth aggregate decreases, the radius of the aggre- 
gate increases, and hence the fractal dimension decreases. 
Figure || shows a typical RA growth cluster at the time 
t = 80tt timesteps. Note that the structure does resem- 
ble meandering vines. Also the radius of the cluster is 
comparable to half of the lattice size (r s=s L/2). 

Figure || is a plot of the average fractal dimension for 
RA clusters as a function of time into the simulation, 
for all three system sizes. The measurements reported 
below are averages over 5 independent realizations for 
the L = 512 system, 10 independent realizations for the 
L = 256 system, and 10 for the L = 128 system (i.e., av- 
erages over either 5 or 10 independently generated, large 
clusters). The data points and standard errors shown in 
the plots are the average and standard error over the set 



of independent realizations. 



The fractal dimension is initially very close to the frac- 
tal dimension for DLA. We then observe a slight increase 
in the fractal dimension as the cluster center begins to an- 
neal (an example is the RA cluster shown in Fig. ||) , then 
a gradual decrease in the fractal dimension until it con- 
verges upon an equilibrium value. The solid line is drawn 
to denote the equilibrium value upon which results for the 
three system sizes are converging, df = 1.6. Using Flory- 
type scaling arguments it has previously been shown that 
a quenched branched polymer obeys the scaling relation- 
ship N ~ R dQ , with d Q = [2(D + 2)]/5 §|. Here R rep- 
resents the characteristic end-to-end distance of a poly- 
mer, and D the dimension of the space. R can be taken 
in direct analogy to I in Eq. ( |l9| ) , which defines the end- 
to-end distance of the window of interest. For D = 2 the 
exponent dQ = 1.6. We should note that an exact result 
was obtained for quenched polymers in D = 2, dQ = 1.64 
pfjj, which is slightly larger. Flory-type scaling has also 
been studied for annealed branched polymers and the 
scaling exponent determined to be dA = (3D + 4)/7 [ j37| . 
For D = 2, dA = 1.43. One might expect to observe 
a crossover from quenched to annealed behavior for the 
equilibrium RA growth clusters as we go from the large 
to the small system sizes, but we did not see this for the 
system sizes investigated. 



Note that the time axis in Fig. ^| is rescaled by L 18 in 
order to match that of Fig. ||. Neither the fractal dimen- 
sion nor the equilibrium temperature exhibit finite size 
effects as far as we can determine within the precision of 
our measurements. 




FIG. 5. A growth cluster grown via the RA model, pic- 
tured at time t = 80tt. The fractal dimension for this cluster 
(df — 1.63±0.02) has seemingly reached the asymptotic value. 
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FIG. 6. The average fractal dimension of the RA growth 
clusters as a function of time into the simulation, for all three 
system sizes. 



VI. DISCUSSION 

We have presented a microscopically reversible model 
which exhibits macroscopic pattern formation. In this 
model, entropy initially grows rapidly with time, and 
then subsequently grows exceedingly slowly — the slow re- 
laxation can be characterized as a quasistatic isothermal 
process. The morphology of the aggregate formed by this 
dynamics changes markedly with time, evolving from a 
pattern having a conformation and fractal dimension sim- 
ilar to that of an irreversible DLA system, to a pattern 
characteristic of a branched quenched random polymer. 

The RA model is an extension of the standard DLA 
model: we model the latent heat released when a gas 
particle aggregates in addition to modeling the gas and 



crystal particles. Since the dynamics is local and mi- 
croscopically reversible, we are realistically modeling the 
flow of heat and the creation of entropy in this system 
and thus we model the thermodynamic behavior of grow- 
ing clusters. 

The model presented is simple and amenable to theo- 
retical analysis. Given the constraints that we set out, 
an even simpler model could be constructed with just a 
single heat bath particle at each site, and a single gas 
particle at each site. In this case, diffusion would be 
performed by block partitioning This simpler model 
would have two fewer bits of persistent state at each lat- 
tice site than the RA model, and would be slightly easier 
to analyze theoretically. It would, however, be less com- 
putationally efficient: for a given lattice size, the volume 
available to the gas and heat particles would be reduced, 
but the computational effort required by each step of the 
simulation would be unchanged. 

There are some simple variants of the RA model which 
might merit study. For example, we have only investi- 
gated situations in which the temperature is set by the 
size of the final aggregate. It would be interesting to 
study morphology in situations where there is indepen- 
dent control of temperature and aggregate size. It would 
also be interesting from a thermodynamic perspective to 
modify the model by introducing a gas-crystal exclusion: 
gas particles would collide with the aggregate, but not 
diffuse over it. Thus there would be an excluded volume 
for the gas particles, and the crystal would do work on 
the gas as it grew. 
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